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Abstract. As an alternative to solving of Poisson equation in Particle-in- 
Cell methods, a new construction of current density exactly satisfying 
continuity equation in finite differences is developed. This procedure 
called density decomposition is proved to be the only possible linear 
procedure for defining the current density associated with the motion of 
a particle. Density decomposition is valid at least for any n-dimensional 
form-factor which is the product of one- dimensional form-factors. The 
algorithm is demonstrated for parabolic spline form-factor. 

1 Introduction 

In the present paper we develope a new procedure called density decomposition for 
obtaining the current density automatically satisfying the continuity equation. 

In the set of Maxwell equations along with hyperbolic equations of wave propa- 
gation we have an equation of elliptic type — Gauss's law, that in terms of electric 
potential (p can be expressed as Poisson equation. In practice Poisson equation is 
used for correction of "potential" part of electric field. 

It is well known that Particle-in-Cell (PIC) method in plasma simulations can be 
implemented without solving Poisson equation for electric field correction. Instead, 
we need the continuity equation (or charge conservation law) in finite differences to 
be satisfyed. 

There are a few methods for satisfying the continuity equation locally — for 
charge and current density associated with each particle, Ref. 0, |]. For this 
purpose authors use special definition for the current density wich is naturally con- 
nected with the change of charge density due to particle motion. Unfortunately, 
these methods are implemented only for simple shapes of particles, for the zero- and 
the first-order form-factors. We present the generalization of these methods, valid 
for a big class of form-factors. Also we have proved that the density decomposition 
is the only possible linear procedure for defining the current density associated with 
the motion of a particle. 

There are another methods for incorporating Gauss's law into Maxwell solver 
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using usual definition of local current density, see ||. 

Very detailed study of PIC method can be found in @, 0, § • The new construc- 
tion will be usefull firstly for overdensed plasma simulation with the paradigm of 
'Clouds-in-Ceir [|. 



2 Continuity equation in finite differences 

Let us consider the local Maxwell solver, wich is equivalent to Finite Difference Time 



Domain (FDTD) method flO 
E"+i - E 



X 
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combined with the particle mover 
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Equations Eqs.(|T]-|) are discreetized Maxwell equations and Eqs.(^^ are leap- 
frog scheme for solving of Newton-Lorentz equations. Here we use dimensionless 
variables defined by transformations t — * 2iTujQ^t, x Aqx, (E, B) — (meCCo'o/e)(E, B), 
where me, e — electron mass and charge, c — speed of light, ujq and Aq — some 
characteristic frequency and length (e.g. the frequency and wavelength of incident 
EM radiation). Index n denotes integer time step and a stands for the number of a 
particle; dt, dx, dy, dz — discreetization of time and space coordinates. 

Different components of electromagnetic fields and charge density p and current 
density J' are defined on different grids, 

E = (-^j,j+l/2,fc+l/2! -^i+l/2,i,fe+l/25 -^i>l/2,j+l/2,A:)5 ^ = (-Bj^;^/2,j,fc5 -^j,j+l/2,fc5 -^i\i,A;+l/2) 5 
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P — Pi+l/2j+l/2,fc+l/2, J — (^iJ+l/2,fc+l/2 5 ^j+l/2J,fc+l/2) ^+1/2 j+l/2,fc) ) 

where z, j, k are integers. Discreet operators in Eqs.([l|-^) are vectors, 

fi+l,j,k fi,j,k fi,j+l,k fi,j,k fi,j,k+l fi,j,k \ 



k 



' dy ^ dz 



V7— -f / f^J'f' fi—i,j,k fi,j,k fi,j—l,k fi,j,k fi,j,k—l \ 

^ = [ Tx ' dy ' d-z J • 

These operators have the next convenient properties 

V- X V+ = V+ X V~ = 0, V" ■ V+ = V+ ■ V" = A, (10) 

where A is discreet Poisson operator in central differences, 

A f. fi—l,j,k '^fi,j,k ~l~ fi+l,j,k fi.j — l,k '^fij.k ~\~ fi,j+l,k fi,j,k—l '^fi,j,k fi.j,k+l 

~ d^^ + di^ + d? 

(11) 

Acting on the Eq.(|l|) by (V~x) and on the Eq.(||) by (V+x), we obtain 

n+l _ n 

+V-- J"+^/^ = 0, (12) 

s = »■ 

It means that if the continuity equation Eq. ([I^ ) is fulfilled then the divergence of E 
is always equal to charge density (Gauss's law), and if the initial discreet divergence 
of B is zero then it remains zero forever. 

Thus, for solving Maxwell equations we need Eqs.(|l]-0) and Eq.(|l2|) with initial 
conditions 

V" ■ E = p and V+ ■ B = at t = . (14) 

Let us consider the continuity equation (or charge conservation law) in finite 
differences 

ri+l/2 _ „ rr\ _ q\ 

A^i+l/2j+l/2,fe+l/2 A^j+l/2j+l/2,fc+l/2 ^ ^ i,j+l/2,k+l/2 <-^i-l,j+l/2,fe+l/2 _^ 

dt dx 

q2 _ '72 '73 _ 7-3 

•-^j+l/2,j,fc+l/2 •-^i+l/2j-l,fc+l/2 •-^j+l/2j+l/2,fc •-^j+l/2J+l/2,fc-l _ p. 

Jy + Jz ^^^^ 
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Further we will drop indices and modificators like ±1/2, where it can not lead to 
an ambiguity. The charge density p is constructed from form-factors of separate 
particles 

Pi,j,k ^ ^ Sij ki^XQ., X/q., ^q.), (16) 
a 

where S is the form-factor (or density) of a particle, 

Xi,Yj,Zk denote coordinates of the grid, {xa,ya,Za) is the location of the particle 
with number a. Here form-factor can be interpreted as a charge density of a single 
particle. So the particle is considered as it would be a charged cloud. Form-factor 
must obey the rule of conservation of full charge which leads to 

SijA^a, Va, Za) = 1, (18) 

where the sum is taken over all grid nodes. 



3 Density decomposition 

Due to linearity of charge conservation law Eq. (|l5|) , it is sufficient to construct 
current density associated with motion of a single particle. 

Let us consider a single particle with form-factor Eq.([T7|) and coordinates (x, y, z). 
We introduce vector W as finite differences of the current density associated with 
particle motion: 



dt ''''''' 

i,j-l,k i,j,k ' 



j2 _ q-2 f.Ty2 

^i.j.k '-^j,7' — l,fc 7i i 



3 3 '^^3 

Then according to charge conservation law, we can write dropping grid indices, 
+ W'^ + = S{x + Ax,y + Ay, z + Az) - S{x, y, z). (20) 
Here {Ax, Ay, Az) is 3-dimensional shift of the particle due to motion. 
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Shift of the particle generates eight functions 

S{x, y,z), S{x + Ax, y, z),S{x,y + Ay, z),S{x, y,z + Az), 
S{x + Ax, y + Ay, z), S{x + Ax, y,z + Az), S{x, y + Ay, z + Az), 

S{x + Ax,y + Ay,z + Az) . (21) 

We will assume that vector W and corresponding current density linearly depends 
from these functions. The base for this assumption is the following. (1) We can 
consider the form-factor as charge density of the particle. If form-factor amplitude 
is increasing, the current density associated with a shift of the form-factor must 
increase proportionally. (2) We can decompose any three-dimensional shift of form- 
factor S{x,y,z) into three one-dimensional shifts: 

S{x + Ax, y + Ay, z + Az) — S{x, y, z) = 
S{x + Ax, y, z) — S{x, y, z) + 
S{x + Ax, y + Ay, z) — S{x + Ax, y, z) + 
^(x + Ax,y + Ay, z + Az) - S{x + Ax,y + Ay, z). (22) 

Currents corresponding to each one- dimensional shift must be additive. 

Let us formulate some conditions directly going form the nature of vector W . 

1. Vector Wlj j^, W^j^i^, Wfj^. is a decomposition of finite difference S'jj,fc(x+Ax, y+ 
Ay, z + Az) - Sij^k{x, y, z), Eq.(||). 

2. If some of shifts Ax, Ay, Az iz zero, the corresponding component W is also 
zero: 

Ax = ^ = 0, Ay = ^ = 0, Az = ^ = 0. 

3. If S{x,y,z) is symmetrical with respect to permutation of {x,y), S{x,y,z) = 
S{y,x,z) and Ax = Ay, then = W"^. The same property is assumed for 
symmetries with respect to permutations of pairs {x,z) and {y,z). 
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Suggestion. There is only one linear combination of eight functions Eg. each 
satisfying Eq.lfTR), that is consistent with properties 1-3; 



+ 
+ 
+ 

+ 

+ 
+ 

+ 
+ 
+ 



= ^S{x + Ax,y + Ay, z + Az) - ]-S{x, y + Ay,z + Az 

+ ^S(x + Ax, y,z + Az) - ^S(x, y,z + Az 
o 

+-S{x + Ax,y + Ay, z) - -S{x, y + Ay,z 

o 

+ ^S{x + Ax, y, z) - ^S{x, y 

W"^ = -S{x + Ax, y + Ay,z + Az) - -S{x + Ax, y,z + Az 
3 3 

1 1 

+ -S{x,y + Ay,z + Az) - -S{x,y,z + Az 

D O 

+-S{x + Ax,y + Ay, z) - -S{x + Ax, y, z 
o o 

+^S{x,y + Ay,z) - -S{x,y 
= ^S{x + Ax,y + Ay, z + Az) - ^S{x + Ax,y + Ay, z 
+^S{x, y + Ay,z + Az) - ^S^x, y + Ay,z 



+-S{x + Ax, y,z + Az) S{x + Ax, y, z 

6 6 



+ ^S{x,y,z + Az) -^S{x,y,z) (23) 



Proof. (Scenario). We can write all the properties 1-3 in the form of linear equa- 



tions with unknown coefficients of eight functions. Remembering Eq.(|T8D we can 
obtain additional equations on coefficients taking sum over all grid points {i,j,k) 
from each linear combination for W. Solving 10 linear equations for all S, we will 
find all the coefficients. Of course, not all eight values Eq.(pTD are independent. 
We have six independend variables x, y, z. Ax, Ay, Az, so in the most general case 
only six values S can be also independend, for example, excluding S{x, y, z) and 
S{x + Ax,y + Ay,z + Az). Among all possible solutions we must left only one, 
which doesn't assume special numerical values for excluded functions. □ 

Taking into account boundary conditions for the current of one particle (van- 
ishing of the current density at nodes far from the form-factor domain), and using 
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Eq.(]T8|) we obtain: 

E = , 

i 

Ewt„fc = o, 

j 

E = • (24) 

k 

Two systems Eq.(^) and Eg . (|2^) define the density decomposition. Solving 
Eq.([19|) with natural boundary condition we obtain the current density associated 
with a single particle motion. 

The condition Eq.(^) can be easily satisfyed if form-factor have a property 
of inheritance in decreasing of the dimension, i.e. if sum of form-factor over any 
dimension is again form-factor but of lower dimension. Formally, it means 

5lf (x,y) = E^S?(^,2/,^), (25) 

k 

where S^^^^ doesn't depend on z and obeys Eg . (pISD automatically. 

There is a big and widely used in PIC codes class of form-factors that have a 
property of inheritance: all form-factors that are the products of one-dimensional 
form-factors, 

SfU^,y,z) = Sl^ix)Sj^iy)Sl^iz). (26) 

Here we use the same symbol for (probably) different one-dimensional form-factors, 
each of them must satisfy conservation of full charge, Eg.(|l^). 

It can be easily proved that density decomposition Eg.(p3D along with Eq.(|2^) 
is the generalization of techniques proposed in 0, ||, |^ . 

4 Computing of the current with second-order 
polynomial form-factor 

In this section we present an algorithm for density decomposition in the case of 
second-order piecewise-polynomial form-factor and discuss a problem of dimension 
reduction. 
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Let us consider well-known one-dimensional form-factor 

S[^\x) = \{^-T{X,-x)) , |X,-x|<l/2, (27) 

which is the second-order spline. The particle is bell-shaped. The correspondent 
3-dimensional form-factor is Eq. (|26|) . 

Now we can formulate a scenario for computing the current density based on den- 
sity decomposition Eg. (pSf) . Suppose we consider a code that uses Finite Difference 
Time Domain (FDTD) technique |T^, where electromagnetic fields and current den- 
sity are defined on different regular grids. Here we do not pretend to show optimized 
or fastest algorithm. 

1. Prepare 15-component array SO containing one-dimensional form-factors cor- 
responding to particle coordinates (xO, yO, zO) with respect to the grid of the 
charge density p: 

SO(^,l) = 5r^(xO),^ = -2,2, 
S0(j,2) = 5f^)(y0),j = -2,2, 

SO(fc, 3) = S'i"^ (zO) , A; = -2, 2 , (28) 

Really, components S0(— 2,m) and S0(2,m) are zero, but we need these addi- 
tional components for further calculations. 

The actual 3-dimensional form-factor is 27-component array 

5^3^) (i, J, k) = SO(z, 1) * SO(j, 2) * SO(A;, 3) . (29) 

2. Using SO or precomputed S*^^^-*, compute the force acting on the particle. Here 
we can use fields spatially averaged to the grid of p or compute additional 
form-factors for each type of grid. Advance particle and compute new particle 
coordinates (xl,yl,zl). Note here that particle shift in any direction must be 
smaller or equal than grid step in this direction, 

xl - xO < dx, yl - yO < dy, zl - zO < dz. (30) 
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3. Using new particle coordinates compute a new array Si containing new form- 
factors: 



Sl(*,l) = 5r^(xl),^ = -2,2, 
Sl(j,2) = ^f^)(yl),j = -2,2, 
Sl(A;,3) = 4'^)(zl),fc = -2,2. 



(31) 



Components Sl(— 2,m) and Sl(— 2,m) are not zero in general, because of 
particle motion. If conditions Eq.(|30D are satisfyed, the array Sl(z, m) doesn't 
have non-zero components out of z = —2, 2. 

4. Compute auxiliary array of differences of new and old form-factors: 



It is possible to use Si for storage of differences. 

5. Compute 125*3-component array containing density decomposition W(i, j, k, m), 
in accordance with Eq.(^). We need so many componets because we have 
current whose components are defined on different regular grids (in FDTD 
technique). 



DS(i,l) = Sl(i,l)-SO(i,l),i 
DS(j, 2) = Sl(j, 2)-S0(J,2),J 
DS(A;,3) = S1(A;,3) -S0(A;,3),A; 



-2,2 



-2,2. 



2,2 



(32) 




(33) 



Of course, this computation is easy to optimize. 
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6. Compute three components of the current density J^^J"^, associated with 
motion of the particle, using Eg. ([T9|) and boundary condition (there is no 
current in nodes far from particle location), 



7^ - 

'•^ i,j,k 


7^ - 

'-'i~l,j,k 


-Qfw(. 


3, k, 1) , 


-72 _ 


't2 _ 
'~'i,j~l,k 




J,fc,2), 


J3 _ 


■ 7^ - 

'^i,j,k-l 


-Q^W(2 

dt ^ 


) J) k, 3) , 



where Q is the charge of the particle. 

7. Add computed contribution from the single particle to array of the current 
density. 

As this algorithm uses only simple polynomes, its accuracy is equivalent to the accu- 
racy of the last digit of numerical representation (e.g. 10~* in SINGLE PRECISION 
4-BYTE data or IQ-^^ in DOUBLE PRECISION 8-BYTE data). 

Suppose we have two-dimensional problem, when all the variables depend on 
{x,y) only. In this case density decomposition Eg . (p3D provides only two first com- 
ponents of the current density. How to construct the third one, in consistency with 
the rest? The simplest idea is to derive the third component from 3-dimensional 
case by reducing the dimension. We can imagine chaines of infinite number of par- 
ticles along 2;-axise. Being projected into (x, ?/)-plane these chaines produces A^ 
2-dimensional particles. Then we can do averaging over z-axise. As a result we will 
obtain first two components of the current density in accordance with Eq.(^), and 
the third component. 

In the particular case of the above algorithm we must change formulae of items 
5 and 6 in the following way: 

W(^, J, 1) = DS(2, 1) * (SO(j, 2) + ^ * DS(j, 2)) , 
W(z, J, 2) = DS(j, 2) * (SO(z, 1) + i * DS(2, 1)) , 
W(^, J, 3) = S0(^, 1) * SO(j, 2) + i * DS(«, 1) * SO(j, 2) + 
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+^ * SO(z, 1) * DS(j, 2) + ^ * DS(z, 1) * DS(j, 2) 



(35) 



^m.-^i = -Qfw(z,j,l), 

^M+i-^M = -Q§W(z,j,2), 

J-,^^. = -QV,W(z,j,3), (36) 

where is the third component of particle velocity. 

As one can see these formulae have an obvious connection with 3D-case, Eqs. 



34) 



5 Conclusion 

In this paper we have developed a construction for a current density, which exactly 
satisfy the charge conservation law and is valid for a wide class of form-factors. It is 
shown that this construction is the only allowed by very natural conditions derived 
from the properties of the current density. An algorithm in the case of second-order 
polynomial form-factor is presented. One can see that this method is not restricted 
by special Maxwell solver, but uses only discreetized continuity equation. These 
teqnique was implemented by author and D.V.Sokolov in three-dimensional and 
two-dimensional PIC codes. 
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